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Abstract An innovative numerical technique is presented to adjust the inflow 
to a supply chain in order to achieve a desired outflow, reducing the costs of 
inventory, or the goods timing in warehouses. 

The supply chain is modelled by a conservation law for the density of processed 
parts coupled to an ODE for the queue buffer occupancy. The control problem 
is stated as the minimization of a cost functional J measuring the queue size 
and the quadratic difference between the outflow and the expected one. The 
main novelty is the extensive use of generalized tangent vectors to a piecewise 
constant control, which represent time shifts of discontinuity points. 
Such method allows convergence results and error estimates for an Upwind- 
Euler steepest descent algorithm, which is also tested by numerical simulations. 
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1 Introduction 

The mathematical modeling of industrial production, as well as the devel- 
opment of techniques for simulation and optimization purposes, is of great 
interest in order to reduce unwanted phenomena (bottlenecks, dead times at 
queues, and so on). Depending on the scale, one can distinguish different mod- 
elling approaches, for instance discrete (Discrete Event Simulations, |14j ) or 
continuous (Differential Equations, [1,2)3,21,22 ). The latter class includes 
models based on partial differential equations ( [51 fTTin2l[TB] ) . For a recent re- 
view see [TO] . 

In this paper, we focus the attention on the continuous model for supply 
chains proposed by Gottlich, Herty and Klar in [16], briefly GHK model. A 
supply chain consists of processors with constant processing rate and a queue 
in front of each processor. The dynamics of parts on a processor is described 
by a conservation law, while the evolution of the queue buffer occupancy is 
given by an ordinary differential equation. The latter is simply determined by 
the difference of fluxes between the preceding and following processors. The 
complete model consists of a coupled PDE-ODE system. 
Various optimal problems, corresponding to different types of controls, have 
been analysed for the GHK model (see [TWrTl[T8lll9[|2"3"] V Typically one may 
consider the input flow to the whole supply chain as a control as well as 
production rates and distribution coefficients in case of supply networks. These 
papers provide a number of results and, in particular, numerical algorithms to 
find the optimal control. In |23j two discretization techniques are compared: 
one consisting in first writing the adjoint system and then discretizing, while 
the second consists in inverting the order, that is first discretizing and then 
writing the adjoint system. All these methods compete with Discrete Events 
ones for numerical accuracy and computational times, see also [TO] . 

In this paper we focus on the optimal control problem, where the control 
is given by the input flow to the supply chain and the cost functional J is the 
sum of time-integral of queues and quadratic distance from a preassigned de- 
sired outflow. In |13j , piecewise constant controls are considered together with 
generalized tangent vectors, which represent time shifts of discontinuities of 
the control. The technique of such generalized tangent vectors was extensively 
used for conservation laws, see [7], and for the case of network models, see [T5l 
[20] . The main result of [13] is the existence of optimal controls. 

The aim of this work is to introduce an innovative numerical approach, 
which builds up on the idea of generalized tangent vectors. Let us first explain 
in rough words the core idea of the approach. A good numerical method, in 
theory and practice, for an optimization problem is often based on a suitable 
choice of "perturbations" . Then one can define critical points, according to 
the chosen perturbations, and numerical algorithms stemming from such def- 
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inition. 

We base our method on perturbations of piecewise constant controls, obtained 
by time shifting the discontinuity points. The advantage of this method is 
twofold: on one side generalized tangent vectors are well suited for Wave 
Front Tracking (briefly WFT, see [51IT5]) algorithms, which allow theoretical 
estimates on convergence rate. On the other side, generalized tangent vec- 
tors have a particularly simple evolution in time, which can be conveniently 
adapted to easily implement able methods as Upwind-Euler (briefly UE, see 
[3]). Moreover, in [5] convergence of the UE algorithm was proved by measur- 
ing the distance with WFT solutions by generalized tangent vectors. 
The discretization of the evolution of generalized tangent vectors allows the 
numerical computation of the cost functional gradient. This can be combined 
with classical steepest descent or more advanced Newton methods for the op- 
timization procedure by iterations. 

Now, the results which complete the picture is the following. To estimate the 
gradient of the cost functional we can use alternatively WFT or UE algorithms 
for the evolution of generalized tangent vectors. Even more, we can measure 
the distance among the two, by again using suitably defined tangent vectors. 
This, in turn, allows to provide convergence results and error estimates. 
Finally simulations are performed to show results for the proposed numerical 
algorithm in some case studies. 

The outline of the paper is the following. In Section [31 we describe the 
GHK model and introduce the optimal control problem. The WFT algorithm 
to construct approximate solutions to the model is illustrated in Section |3l 
together with the definition and evolution of generalized tangent vectors. Then 
Section @] describes the UE numerical algorithm for solutions to the coupled 
ODE-PDE system and also the numerics for generalized tangent vectors and 
cost functional derivative. Convergence results and rates for our method are 
then given in Section [5] Finally some numerical tests are discussed in Section 

m 

2 An optimal control problem for supply chains 

A supply chain consists of consecutive suppliers. Each supplier is composed 
of a processor for parts assembling and construction and a queue, located in 
front of the processor, for unprocessed parts. Formally we have the following 
definition. 

Definition 1 A supply chain consists of a finite sequence of consecutive pro- 
cessors Ij, j G J = {1, . . . , P} and queues in front of each processor, except 
the first. Thus the supply chain is given by a graph G = (V, Sf) with arcs rep- 
resenting processors and vertices, in V = {1, . . . , P — 1}, representing queues. 
Each processor is parametrized by a bounded closed interval Ij = [a,j, bj], with 
bj-i = aj,j = 2, ...,P. 

The maximal processing rate /x^, and the processing velocity, given by 
Vj = Lj/Tj with Tj and Lj = bj — aj the processing time and the length of 
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the j-th processor, are constant parameters for each arc. The dynamics of the 
j-th processor is given by a conservation law 



d t pj (x, t) + d x min {pj,Vjpj {x, t)} 

p j (x,0)=p j , (x), Pj(aj,t) 



0, Va; G [dj, bj] , t G 



(1) 



where pj G [0, pj lax ] is the unknown function, representing the density of parts, 
while the initial datum pj$ and the inflow fj t i nc (t) are to be assigned. An input 
profile u(t) on the left boundary {(ai, t) : t G K} is given for the first arc of the 
supply chain. Each queue buffer occupancy is modelled as a time-dependent 
function t — > qj(t), satisfying the following equation: 



= fj-l (pj-i (bj-i,t)) - fj, 



i = 2,...,p, 



(2) 



where the first term is defined by the trace of Pj-i (which is assumed to be of 
bounded variation on the x variable), while the second is defined by: 



min{/j_i (pj-i (bj-i,t)) ,p,j} if qj (t) = 0, 
Pj if qj (t) > 0. 



(3) 



This allows for the following interpretation: If the outgoing buffer is empty, 
we process as many parts as possible but at most pj. If the buffering queue 
contains parts, then we process at the maximal possible rate, namely again 
Pj . Finally, the supply chain model is a coupled system of partial and ordinary 
differential equations given by 



' d t pj (x, t) + d x min {pj,Vjpj (x, t)} = j = 1, 

Qjit) = fj-l (Pj-i 0>j-i,t)) - fj,inc(t) j = 2, 

qj (0) = q J:0 j = 2, 

Pj (x, 0) = p jfi (x) j = l, 

/9 .( a . )t) = teW j = h 

> /l,tnc(t) = u{t) 

where fj.i nc is given by ([3]), for j = 2, ...,P. 

Fixed a time horizon [0, T], define the cost functional: 



■;P, 
■;P, 
..,P, 

■;P, 
..,P, 



(4) 



J(u)=Y] [ ai(t) qj {t)dt+ [ a 2 (t) [v P ■ p P (b P ,t)) - ^{t)\ 2 dt 

„-_n JO JO 



J'=2 

Ji(u) + J 2 (u), 



where a x G ^((O.T), [0,+oo)) 5 a 2 G Lip((0,T), [0,+oo)) (the space of Lips- 
chitz continuous functions) are weight functions, (/3j, qj) is the solution to Q 
for the control u, vp ■ pp(bp,t)) represents the outflow of the supply chain 
(assuming the density level is below pp), while ip(t) G Lip((0,T), [0, +oo)) is 
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a pre-assigned desired outflow. Given C > 0, we consider the minimization 
problem 

min J (u) (6) 

uGUc 

where U c = {u : [0,T] -> [0,jUi] ; « measurable, T.V.(u) < C*} (with T.V. in- 
dicating the total variation). In other words, we want to minimize the queues 
length and the distance between the exiting flow and the pre-assigned flow 
ip{t), using the supply chain input u as control. 

In [T3] , the existence of an optimal control was proved for a general prob- 
lem, which includes the case (JU)-©: 

Theorem 1 (see flSjj) Consider the optimal control problem O), (0). If J is 
lower semicontinuous for the L 1 norm, then there exists an optimal control. 

Our aim is now to provide a new approach to solve (U)-© numerically. 
The key idea is to focus on piecewise constant controls and perturb the posi- 
tion of discontinuity points. The procedure corresponds to define (generalized) 
tangent vectors to u (in the spirit of [7]). We can then take advantage of the 
knowledge of time evolution of such tangent vectors, developed in [20] . 
More precisely, we start giving the following: 

Definition 2 Wc indicate by U C Uc the set of Piecewise Constant controls. 
For every «£l/wc indicate by = Tfc(it), k = 1, . . . , S(u), the discontinuity 
points of u. 

We are now ready to define the perturbation to a piecewise constant con- 
trol: 

Definition 3 Given u eW, a tangent vector to u is a vector £ = (£i, . . . , £s(u)) £ 
U<5(«) representing shifts of discontinuities. The norm of the tangent vector is 
defined as: 

5(u) 

IKII = El&l>fa+)-«fa-)|. 

Assume now for simplicity that t\ > 0, Tgr u ) < T, and set tq — 0, £o = 0, 
T S(u)+i — T, £,s(u)+i = 0. Then given a tangent vector £ to u, for every e 
sufficiently small we define the infinitesimal displacement as: 

S(u) 

u e = E X[T h +s£ k ,T k+1 +e( h+ i[U(T-k+), (7) 
k=0 

where x is the indicator function. In other words u e is obtained from u by 
shifting the discontinuity points of e£ (see Figured]). 

Remark 1 The notion of tangent vector to piecewise constant function was in- 
troduced in [5] and used in [7] to prove uniqueness and continuous dependence 
of solutions to systems of conservation laws. 
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r, r 2 r 3 

Fig. 1 Shifting of the discontinuity points of the control u. 

More precisely, one defines a distance of Finsler type among piecewise con- 
stant functions by considering paths which admit tangent vectors. Then, by 
density and using the usual L 1 metric one can extend the metric to the whole 
L . Finally one may study the evolution of tangent vectors and, in particular, 
estimates on their norms. 

This same technique was generalized to flows on networks, see |15) . and to the 
GHK supply chain model, see [20] . 

In next sections we are going to define numerical schemes for the solution 
of and for the evolution of tangent vectors. The latter, in turn, will provide 
the information for the computation of numerical gradient of the cost func- 
tional J. 

The evolution of tangent vectors is particularly clear for the theoretical nu- 
merical scheme given by the WFT algorithm. This was used in [5] to prove 
convergence of an Upwind-Euler scheme and in [T3] to obtain the existence of 
an optimal control for ((U)-®. We thus recall briefly the WFT algorithm and 
the evolution of tangent vectors along approximate solutions constructed via 
the WFT algorithm. 



3 The Wave Front Tracking algorithm 

In this section we explain how to construct piecewise constant approximate 
solutions to dl]) by WFT method, see [B] for details. 

Given a discretization parameter a and initial conditions pj t o in BV, the 
space of bounded variation functions, a WFT approximate solution is con- 
structed by a procedure sketched by the following steps: 

— Approximate the initial datum by a piecewise constant function (with dis- 
cretization parameter a) and solve the Riemann Problems (RPs) corre- 
sponding to discontinuities of the approximation. In RPs solutions approx- 
imate rarefactions by rarefaction shocks of size a; 
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— Use the piecewise constant solution obtained piecing together the solutions 
to RPs up to the first time of interaction of two shocks; 

— Then solve the new RP created by interaction of waves and prolong the 
solution up to next interaction time, and so on. 

To ensure the feasibility of such construction and the convergence of WFT 
approximate solutions to a weak solution as a — >■ 0, it is enough to control 
the number of waves and interactions and the BV norm. This is easily done in 
scalar case since both the number of waves and the BV norm are decreasing 
in time (see [B] for details). 

For our system, we need also to approximate the queue evolution. For this we 
compute the exact solutions to (|2|) (see [20] for BV estimates for the complete 
PDE-ODE model Q). 

Notice that, as soon as a boundary datum will achieve a value below p,j, 
then in finite time all values above pj will disappear from the j-th processor, 
see also [3U]- Therefore, for simplicity, we will assume 

(HI) Pjfi{x) < fj,j for all j E J. 

Then the same inequality will be satisfied for all times. In this case solutions 
to RPs are particularly simple, indeed the conservation law is linear, thus given 
some Riemann data (/?_, p + ) on the j-th processor, the solution is always given 
by a shock travelling with velocity Vj . 

3.1 Tangent vectors evolution 

The infinitesimal displacement of each discontinuity of the control u produces 
changes in the whole supply chain, whose effects are visible both on processors 
and on queues. In fact, every shift £ generates shifts on the densities and shifts 
on the queues, which spread along the whole supply chain. 
Since the control u is piecewise constant, the solution (pj,qj) to ((4]) is such 
that pj is piecewise constant and qj is piecewise linear. A tangent vector to 
the solution (pj,qj) is given by: 



where j3 runs over the set of discontinuities of pj, are the shifts of the 
discontinuities, while r/j is the shift of the queue buffer occupancy qj. The 
norm of a tangent vector is given by: 



Notice that this is compatible with the definition of norm of tangent vector 




(8) 




8 



Ciro D'Apice et al. 



to a control. Because of assumption (HI), we have no wave interaction inside 
the processors. Therefore, densities and queues shifts remain constant for al- 
most all times and change only at those times at which one of the following 
interactions occurs: 

a) interaction of a density wave with a queue; 

b) emptying of the queue. 

We now provide formulas for such changes. Assume a wave with shift 
interact with the j-th queue and let t be the interaction time. We use the 
letters + and — to indicate quantities before and after t, respectively. So, we 
indicate with pj and p^ the densities on the processor Ij before and after an 
interaction occurs and similarly for Ij-i- Also we use 13 £j to denote the shift on 
the processor Ij and with @ rj~ and ^rf^ the shifts on the queue qj, respectively 
before and after the interaction. Consider the case a) and distinguish two 
subcases: 

a.l) qj (i) = 0; 
a.2) qj(t) > 0. 

In case a.l) we have to further distinguish two subcases: 
a.1.1) if < "j-iPj-i < A*j> thcn % = TT^T^Ci-i and f3 r]j = 

o = e v + ; 

a.1.2) if pj-xpjL! > ^,then%- = , and^+ = + 

%■ 

In case a.2) we have: ^ = 0, 13 = ^•_i(pj"_ 1 - Pj_ t ) + f3 Vj- 

R — 

Finally in case b) we get: = 0, = and ^£ 



Using the above notations, indicate by ^£p the shift to a generic discon- 
tinuity of pp and by ^Pp, respectively 13 p~p~, the value of pp on the right, 
respectively left, of the discontinuity. The following holds: 

Proposition 1 Consider a control u G hi and a tangent vector £ £ R 5 (") to 
u. The gradient of the cost functional J with respect to £ is given by: 

V C J(«)= V a 1 (t) Vj (t)dt + Va 2 (t> P (^+ !3 t P A{' 3 pp) 
j Jo 

= Y™ FT + F 2 ^ T , (9) 

w/iere /A^pp) = 13 pp — 13 Pp andt 13 is the interaction time of the discontinuity 
indexed by (3 with bp, the right extreme of the supply chain. 

Proof We have J(u) = lim e ^ ^Xl^z^M By definition of the functional J, 
the infinitesimal change Y^ FT in J x due to the infinitesimal displacement e£ 
of the control u is 



Y WFT = £ ^2 f ai (t%j(t)dt, (10) 

j Jo 
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while the infinitesimal change Y^ FT in J2 is 

yWFT =£ J2a 2 (tP) V P [(^+) 2 - ^ Pp ) 2 (*p+ - *pp) 



thus we conclude. 



(11) 



4 Steepest descent for the Upwind-Euler scheme 

In this section we introduce first an Upwind-Euler scheme for the system 
(UJ and then a numerical scheme for the evolution of the tangent vectors 
to a solution to the PDE-ODE model. From the latter we will be able to 
compute numerically the derivative of the cost functional with respect to the 
discontinuities of the input flow. This, in turn, will be used in steepest descent 
methods to find the optimal control. 

For a general introduction to numerical schemes for conservation laws we refer 
to [24 , while for optimization algorithms to [I]. 
For simplicity we assume: 

(H2) The lengths Lj are rationally dependent. 

Assumption (H2) allows us to use a unique space mesh for all processors 
Ij , j = 1, . . . , P. Indeed there exists A so that all Lj are multiple of A and we 
will always use time and space meshes dividing A. 

Remark 2 It is possible to choose different space and/or time grid meshes for 
different processors. This is necessary in the general case in which the lengths 
of arcs have not rational ratios. Details can be found in [S]. 

In next section we report briefly the Upwind-Euler method, analysed in [5] 
to construct numerical solutions to the supply chain model 



4.1 Upwind-Euler scheme for supply chains 

Given a space mesh Ax, for each processor lj, we set Atj = Ax/vj and define 
a numerical grid of [0, Lj] x [0, T] by: 

— (xi,t n )j — (iAx, nAtj), i — 0, Nj, n — 0, Mj are the grid points; 

— 3 pf is the value taken by the approximated density at the point {xi,t n )p 

— is the value taken by the approximate queue buffer occupancy at time 
t n . 

The Upwind method reads: 
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where j e J , i = 0, Nj and n = 0, My. Notice that the CFL condition is 
given by Atj < — , and thus holds true. The explicit Euler method is given 
by: 

q] +1 = q] + Atj (/;_! - ff Mc ) , j e J\ {1} , n = 0, M h (13) 
where f™_ 1 needs to be defined and 

Jj,inc \ k tj / 1 \ rv » ' 

I Mj ?" (*) > 0. 

Now, if Atj < Atj we set: 

M(n)—m(n) — l 7 

;=i i=i 

(15) 

where m (n) and M (n) are defined as: 

m (n) = sup {to : mAtj-i < nAtj} , 

M (n) = inf {M : MAtj^ > (n + 1) . 
Otherwise, that is if Atj-\ > Atj, we set: 

/r, /, W 1 ^)- as) 

where ['J indicates the floor function. Boundary data are treated using ghost 
cells and the expression of inflows given by (|14p . The convergence of the scheme 
has been proved in [5] using a comparison with WFT approximate solutions. 



4.2 Numerics for tangent vectors and cost functional 

We first completely discretize the control space via the time mesh At: 

U A t = {u GW : T k (u) = n(u, k) At, n(k, u) e N, k = 1, . . . , S(u)}. (17) 

Now for every u € U A t we consider shifts £ so that the obtained time-shifted 
control is still in hi At- Then every is necessarily a multiple of At. Hence 
from now on we will restrict to the case: 

ik = ±At, k = l,...,6(u). (18) 

For a generic processor Ij and a discontinuity point of the control, 
we denote by " and k ^r\ 71 the approximations of k £j(xi,t n ), and k r]j(t n ), 
respectively. We define such approximations by a recursive procedure explained 
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in the following. 

We initialize the tangent vector approximations by setting: 

= 0, farn=l,...,n(fc-l,u), j = l,...,P, 

fc 'V = 0, j = l,...,P. 

The definition of M^C ' 11 ) reflects the fact that the shift provokes a shift 
of the wave generated on the first processor. 

Now, the evolution of approximations of tangent vectors to p inside processors 
is simply given by: 

fcjVn+l k,j cn 

On the other side, the approximation of £ and 77 influence each other at inter- 
action times with queues. More precisely, we consider the four cases described 
in Section IQl and get: 



.1.1): fe >V +1 = 0, fc ^r +1 = ^ fcj '- 1 ^ J _ 1 ; 

n /c,j' n+1 k 

N4-1' " ~ " ^Nj 



a 1 2V ' c ,.7£' rl +l = Vl k,j — lpn k,j n+l _ 



a.2): ' 1 •). fc ^" +1 = fc ^- 1 a,_ a (;' Vy' - ' Vy : ) + ? ' "<'/": 
b): = 0, = 0, fe >->£" +1 = • 



Notice that this is compatible with the evolution of tangent vectors along 
WFT approximate solutions and also the norm of approximate tangent vectors, 
in the sense of ([5]), is conserved. 

Now we are ready to compute numerical approximations for J. We de- 
note by fe ^Y 1 ™, respectively k Y£ l , the numerical approximations of the A:-th 
component of Y^ FT , respectively Y™ FT , as defined in © on processor Ij at 
time t n . 

We initialize such approximation by setting: 

Wyo = 0j k Y o = Qj j = 1, . . . , F, fc = 1, . . . , S(u). (20) 

Now the evolution is determined by the following simple rules. For Yi, if > 
0, then we set 

fe,jyn+l = ^Y? + ai(t n ) k ' j V n At, 
while for q" +1 =0we distinguish two subcases: 

- if = 0, then = ^jyn. 

- if qf > 0, then k ^Y™ +1 = + 1 ai (jn) fc,j£«+i k^n 
For Y" 2 we set: 



k\rn-\-l k\rn , /.\ / /P n / /j.7l\\ (P n+l / \ k,P <-r, 

Y 2 = Y 2 + a 2 {t)v P ^{ p Np -ijj{t )) -( p^ p -tp(t )) j ■ £ Nf 
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A steepest descent algorithm, denoting whit $ the iteration step, is defined 
by setting 



where hg is a coefficient to be suitably chosen. More precisely the parameter 
hg may be chosen to solve an optimization problem to get specific schemes, 
see @]. 

5 Convergence and error estimates 

In this section we will provide convergence results and error estimates for 
the Upwind-Euler steepest descent scheme illustrated in Section SJ We will 
make use of two natural parameters v 6 N and 6 £ N, the first referring to 
the Upwind-Euler (and WFT) discretization, while the second indicating the 
iterative step of the steepest descent algorithm. 

We fix an initial space mesh Ax (of which A is a multiple, see (H2)) and 
define Ax v = 1~ v Ax. On each processor Ij the time mesh is set as: 



thus granting the CFL condition. Obviously Atj lU — > 0, Ax v — > as v tends 
to +oo. 

The initial datum (see Q) is sampled in the following way: 



A control function u Ul g will be defined by the iteration step of the steepest 
descent algorithm starting from a fixed u Vj q G UAti „ ■ We will denote the 
discontinuity points of u u ,e by t^' 6 , for k = 1, . . . , 8{u u ^g). 
Notice approximations can be constructed in such a way that 8{u u fi) — > +oo 
as is y +oo. 

For simplicity v,e (p : q) UE will denote the numerical solution generated by 
the Upwind-Euler scheme, i.e. the collection ly ' e (- J p", g") for j — 1,...,P, 
n = l,...,M!f,i = l,...,N!f. Similarly "• e (£ : , if) UE will denote the numerical 
tangent vectors computed as in Section H i.e. the collection ufi { k ' 3 Ci, k,3r n n ) 
for j = 1, . . . ,P, n = 1, . . . ,Mjf, i = 1, . . . , Nf, k = 1, . . . , S(u v , e ). 
We will also use the symbols ' (p, q) WFT , respectively v ' r]) WFT , to indi- 
cate the solutions, respectively tangent vectors, produced by the Wave Front 
Tracking algorithm. 

Now define: 



{"' j p n ) = X! X[a j +i2- v Ax j ,a i +(i+i)2- v Ax j l> ( 22 ) 





(21) 



L j /2~' J Ax j -l 
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and 

7T PL (t) = « q *+(t-nAt jiV ){ v q™ +1 - u qj) for i e [n^ >5 (n+l)^>[. 

Notice that npc, respectively ttpl, are operators taking values on the space 
of piecewise continuous, respectively piecewise linear, functions. 
In [5] it was proved the following: 

Theorem 2 Assume that (HI), (H2), fHp hold true and that pj^ are of 
bounded variation. Then ub \p,q) WFT is approximated by u '° (irpc(p),TTpL(q)) UE , 
more precisely: 

r e p wFT {t) _ •>,' {vpc {p UE W)\\» + r e qWF T {t) u,e {7rpL{qU E )m 
<2~" K(T.V.( Pjt0 ), N ), (23) 

where K is a suitable constant depending only on the total variation of pj^ 
and the values fij, j = 1, . . . , P. 

Moreover "' e (irpc(p), ^pl{o)) UE converges to a solution of ^ with a conver- 
gence rate as in 

The main idea behind Theorem [2] is to use tangent vectors to estimate the 
distance among WFT and UE approximate solutions. In the same fashion we 
are now going to estimate the distance among tangent vectors computed via 
WFT and UE schemes. 

In each processor Ij, the WFT wave shifts ^^j{t n ) are approximated by 
the UE shifts k *£? in the following sense. Each shift ^£,j{t n ) corresponds to 
one or more non vanishing values of fc,J £™, whose position is possibly shifted 
by some amount. By splitting the shift ^£,j{t n ) in one or more pieces, whose 
sum is equal to P£j{t n ), we can define the difference with the UE generated 
shifts by tangent vectors ^(j(t n ), which represent the space distance among 
the location of the discontinuity (3 of pj at time t n and the location of the non 
vanishing values of i.e. a,j + i Ax v . 

For instance in the WFT algorithm a shift of the time discontinuity Tfc, 
gives rise to a shift £ = i>i£jt on processor 1\ at time t% of the wave generated by 
the discontinuity Tk ■ Similarly the UE algorithm will generate a non vanishing 
value > u ) — v\^k- Therefore, we have Ci( T fc) — Ci( n (&> u)Ati tl/ ) = 0, 

indeed the WFT and UE shifts coincide. 

Analougsly we define the tangent vectors ij to the queue shifts rjj generated 
by WFT to recover those generated by UE. 

We define the norm of the tangent vectors (C, i) in the following way: 

ikc.oii = EE f&i \M\ +E w fei- ( 24 ) 

In particular we have: 

||(C,0(0)||=0. (25) 

Such tangent vectors (£, l) evolve using the same rules of (£, rf) for the 
Wave Front Tracking algorithm, see Section 13.11 In particular their norm may 
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increase because of interactions with queues. Moreover, the norm may also 
increase because of errors produced by the use of different time meshes among 
different processors, see formula (|15[l and (|16|) . Summarizing, to estimate the 
increase in (£, l) norm we have to consider: 

i) approximation errors occurring at interaction of waves with queues and at 
time in which a queue empties; 

ii) the approximation errors due to different time discretizations At^ v . 

Let us start by considering the case i) assuming there is no error approx- 
imation because of different time meshes. Since the evolution of (£, i) follows 
the same rules as those of (£, rj) for WFT, we can use the same estimates es- 
tablished in [9]. Referring to formulas (29), (30) and (31) of [9], and using the 
symbols ± to indicate quantities before and after the interaction, we get: 



IKCO+II < ^l^' 1 / ) IKC:0-ll+%h7l^lA°;l, (26) 

where j is the queue involved in the interaction and Apj is the jump in p of 
the wave exiting to processor Ij . More precisely, the first term on the right of 
(f2"6"| is sufficient for the case of interaction of a wave with a queue, see also 
(29) and (30) of [9]. In case of emptying of the queue qj, the first term takes 
into account the evolution of (rj, i), as described in Section f3.ll The second 
term takes into account the additional shifts, provoked by the fact that the 
WFT produces a wave in Ij at a time which may be not a multiple of Atj tV 
(as it happens for all waves of the UE algorithm), see also formula (31) of [9]. 

Regarding ii), we refer to formulas (TT5j) and ([TBI . In case of Atj-i tU > Atj t „, 
then no additional error occurs. Otherwise, the error is estimated by: 

IKCO+II < IKCO-II +v j At j - 1 , v I \Apj-i\, 

where is the shift of the interacting wave and Apj—\ is the jump in p of 
the interacting wave. 

Recalling (|25p . from the above estimates we get the following: 

sup ||(C,0(t)||<IT mflX i — '4 
te[o,T] fJ- 2 [Vj-i J 

■maxvjAt j<v - sup \\ v ' 6 ^,v) WFT {t)\\ • sup T,V,{ V > 6 p WFT {t)). (27) 
J te[o,T] te[o,T] 

Now in [3D] it was proved (Lemma 2.7) that the norm of tangent vectors are 
decreasing along WFT solutions. Moreover, the tangent vectors for the WFT 
solutions satisfy the same initial estimate as (|19D . thus we get: 

sup || vfi {t,rj) WFT {t)\\ < || v ' e (£,r,) WFT (0)\\ - «i Ah,* - Ax v = 2^ Ax. 
te[o,T] 

(28) 
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From (2.10a) of [20J and (EI]), we get: 

sup T.V.(»> e p WFT (t)) < T.V.(»> e p WFT (0)) + \\d t (^ e q WFT 
te[o,T] 

• ^7.l .; /); „: ■ £ //, ,,, , . (29) 

3 3 

Since maxj VjAt^ v = Ax v = 1~ v Ax, we get: 

Theorem 3 Assume that (HI), (H2), 121)) hold true and that pj$ are of 
bounded variation. Then the following estimate holds true: 

sup || (COW II 

te[o,T] 

< J]max|-^-,lJ U^T.V.^o) • £ p J //„ , j . 

We are now ready to prove the following: 

Proposition 2 Assume that (HI), (H2), 121)) hold true and that pjfi are of 
bounded variation. Let Y^ JE , i — 1,2, indicate the numerical approximation 
of Yi via the Upwind-Euler steepest descent scheme of Section [^] Then there 
exists a constant K\ depending only on the data of the problem: 

Ki = Ki (||ai||ii, Lip(a 2 ), \\a 2 \\oo,Lip(rp), ||^||oo,«p, Vj>T.V.(pj,o)) , 

where Lip(-) indicates the Lipschitz constant of a function, \\ ■ ||oo the L°° norm 
and j — 1, . . . , P, such that the following estimates hold: 

\\Y^ FT -npciY^W <K X sup ||(C,0(*)||, 

te[o,T] 

where itpc is the projection over piecewise constant functions as in \22)) and 

\\Y™ FT -YV E \\< Kl sup ||(C, t )(t)||. 

te[o,T] 

Proof From we have: 

^Y?™ -*PC<y? E )\\<\\ax\\v sup \\{Q,C){t)l 



te[o,T] 

and 



lY WFT_ Y UEu < 



[Lip(a 2 )vp(U\\ 00 + 2pp) sup \r 6 (t,v) WFT (t)\\ sup T.V.{ v >°p WFT (t)) 
\ te[o,r] *e[o,T] 

+ \\a 2 \\ 00 vpLip(i>) sup r 9 (i,v) WFT {t)\\ sup T.V.C> e p WFT (t))) 
te[o,T] te[o,T] ) 

■ sup ||(C,0(*)ll- 
fe[o,T] 

By ([28]) and we conclude. 
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We are now ready to state our main convergence result: 

Theorem 4 Assume that (HI), (H2), 121)) hold true and that pj.o are of 
bounded variation. Fixing v, if 5{u v _g) = S , hg > h > 0, for all 8, and 
t^' 9 — ¥ ffc as 8 — >• +oo, then u Vi g strongly converges in L 1 to some u 6 Uc as 
8 — > +oo and 

V e J(u) < K 2 2~ v Ax 

where K 2 depends only on the data of the problem as for K\ (of Proposition 
0) and on Vj , j = 1, . . . , P. 

Proof Clearly u v fi strongly converges in L 1 to some u, moreover by Helly 
theorem u e Uc- 

Now from Proposition [Q we have that V ( J(u„ ld ) = Y^' FT + Y 2 WFT . The 
tangent vectors (£, rj) WFT converge as 8 — > +oo, in the sense that the 
positions and values of shifts converge. Then we get that J{u v ^e) converges 
to V ? J(u). 

Now since t£ 9 — > ?k we have that Y^ E and Y^ E are bounded by At\ iV (times 
a constant depending only on the data of the problem as for K\). Thus we 
conclude by Theorem [3] and Proposition [5] 

6 Simulations 

In this Section we test the Euler-Upwind steepest descent algorithm on two 
test cases. 

Consider first a supply chain characterized by 11 arcs with the following 
input function: 

{ui < t < n, 
U 2 Tl < t < T 2 , 
U 3 T 2 <t<T. 

We assume that the supply chain is initially empty and set Vj = 1 Vj = 1, 11 
and 

/ii = 200, /12 = 75, /13 — 100, /i4 = 65, /15 = 150, 

/i 6 = 75, /Lt 7 = 30, ^ 8 = 100, fig = 80, /iio = 100, /in = 120. 

For simplicity we also set ot\ = 1 and a 2 = and analyze two different cases: 

Case a) u\ — 90, u 2 = 100, 1*3 = 125, T= 10, initial values (ti,t 2 ) ~ (1,3); 
Case b) m = 100, u 2 = 90, u 3 = 125,T= 10, initial values (ti,t 2 ) = (4,5). 

As time and space grid meshes we choose Ax — 0.02 and At — 0.016 so as 
to satisfy the CFL condition. We use the condition that J remains unchanged 
for five runnings of the algorithm as forced stop criterion. The times n and 
t 2 found by the algorithm are: 

Case a) t\ ~ 8.98, t 2 ~ 9.14: as expected both t\ and t 2 run toward T; in 
fact, in order to minimize the queues, the optimization algorithm tends to 
reduce the inflow levels which increase the queues (i. e. u 2 and U3). 
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Case b) n = 0, T2 ~ 9.05: ri runs to zero and Ti runs toward T; as in the 
previous case, the optimization algorithm works to reduce the inflow levels 
which lead to queues increasing (i. e. u\ and U3). 

In Table [T] we report the numerical values of t\ , t-i and J at each iteration of 
the steepest descent algorithm for Case a). 



Table 1 Case a: J versus n , T2 in 42 iterations of the steepest descent algorithm. 



Iteration 


Tl 


T2 


Ji 


Iteration 


7~1 


T2 


Ji 


1 


1 


3 


117799.059 


22 


8.873 


9 


076 


69504.563 


2 


2.595 


5.988 


89474.759 


23 


8.887 


9 


084 


69503.235 


3 


3.870 


7.480 


79616.859 


24 


8.897 


9 


093 


69501.100 


1 


4.893 


8.228 


75291.519 


25 


8.906 


9 


099 


69501.100 


5 


5.710 


8.600 


73053.659 


26 


8.915 


9 


105 


69500.948 


6 


6.365 


8.788 


71737.0579 


27 


8.923 


9 


110 


69500.948 


7 


6.888 


8.880 


70923.172 


28 


8.931 


9 


115 


69500.073 


8 


7.306 


8.926 


70419.523 


29 


8.937 


9 


119 


69500.073 


9 


7.640 


8.957 


70092.102 


30 


8.942 


9 


124 


69500.061 


10 


7.907 


8.977 


69879.287 


31 


8.949 


9 


1276 


69499.327 


11 


8.120 


8.987 


69749.807 


32 


8.954 


9 


130 


69499.327 


12 


8.291 


8.998 


69660.449 


33 


8.959 


9 


134 


69499.327 


13 


8.427 


9.005 


69608.289 


34 


8.964 


9 


137 


69498.722 


11 


8.538 


9.012 


69570.370 


35 


8.968 


9 


139 


69498.722 


15 


8.625 


9.017 


69544.690 


36 


8.972 


9 


140 


69498.722 


16 


8.694 


9.022 


69530.770 


37 


8.976 


8 


972 


69498.245 


17 


8.748 


9.028 


69521.531 


38 


8.978 


9 


144 


69498.245 


18 


8.790 


9.035 


69514.2101 


39 


8.981 


9 


146 


69498.245 


19 


8.818 


9.047 


69509.809 


40 


8.984 


9 


147 


69498.245 


20 


8.839 


9.058 


69507.807 


41 


8.986 


9 


149 


69498.245 


21 


8.858 


9.066 


69506.164 


42 


8.989 


9 


150 


69498.245 



Figure^ respectively [31 depicts the values assumed by J, respectively (ri , 
T2), at each iteration step for Case a). 



120 000 h. 



— 1 — ^ Iteration 

40 



Fig. 



2 Supply chain with 11 arcs, Case a. J versus iteration steps. 
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8 - 



4 - 



2 4 6 8 

Fig. 3 Supply chain with 11 arcs, Case a. Path followed by the steepest descent algorithm 
in the plane (ri , ra) . 



The behaviour of the cost functional J in the plane (t%,T2), is reported for 
Case b) in Figure HI to confirm the goodness of the steepest descent algorithm. 
Notice that since J decreases when the number of iteration increases, the 
updating of t\ and values allows an effective decrement of queues at supply 
chain nodes. 




Fig. 4 Supply chain with 11 arcs, Case b. Behaviour of J\ in the plane (ti,T2). 



We now analyze now a supply chain with 2 arcs, maximal processing rates 
jiii = 200, /i2 = 75 and lengths and processing rates of each processor equal to 
1. We assume that processors and queues arc empty at t — 0, i.e. pj o (x) = 0, 
Vie [0, 1] , j = 1, 2, <7i(0) = 0. The levels of the input flow are u\ = 100, u-i — 
80, U3 = 50. The total simulation time is T = 20 and numerical approximations 
are made with Ax = 0.02, At — 0.016. The aim is to optimize J = J\ + J2,with 
a.\ = ct2 = 0.5, i.e. minimize the queue handling a pre-assigned piecewise 
constant outflow: 
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i> = r/>{t) 



100 < t < 10, 
75 10 < t < T. 



Figures [5] and [5] show the values assumed by J at each iteration step and the 
"path" followed by the steepest descent algorithm in the plane (n, r 2 ), starting 
with initial searching point (ti,t 2 ) = (5,12). We observe that according to 
the aim of minimizing the queue, J is a decreasing function and, moreover, as 
expected the flux on the last arc is equal to the final outflow value. Finally J 
is minimized (its value is zero) by (£1,^2) = (0,0). 



J 

200 000 r 



Iteration 



2 4 6 8 10 12 14 

Fig. 5 Supply chain with 2 arcs: J versus iteration steps. 



T2 

12 - 



Fig. 6 Supply chain with 2 arcs: "path" followed by the steepest descent algorithm in the 
plane (ti,t 2 ). 
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